clearvars -except b

load tprob_data 

tempind = (tprob_t>2);
tprob_m = tprob_m(tempind);
v_crime = v_crime(tempind);
Lv_crime = Lv_crime(tempind);
tprob_t = tprob_t(tempind);
LLv_crime = LLv_crime(tempind);


num_coef = 4;

M = max(tprob_m);
pxest = zeros(M,num_coef);sigrhoest=zeros(M,1);

for m = 1:M
    index = find(tprob_m==m);
    cons = ones(length(index),1);
    [temp ,~,r,~,stats] = regress(v_crime(index), [cons Lv_crime(index) tprob_t(index) LLv_crime(index)]);
    pxest(m,:) = temp;
end
pxest = pxest';


T = 7; 

yvals = [1:1:19]';
y=1;
gamma = ones(M,1);
beta = 0.95;
xl = 2;

v_crime_bar_big = 999*ones(M,2);
for val = 1:2
    x = val-1;
    v_crime_bar = x*ones(M,1);
extp=x*ones(M,1);
extpL=xl*ones(M,1);
for m=1:M
B = 1;
for t=2:7 
extpL(m) =  extp(m);    
extp(m) =  pxest(1,m) + pxest(2,m)*extp(m)  + pxest(3,m)*(yvals(y)+t-1) + pxest(4,m)*extpL(m) ;


temp = beta^(t-1);
temp2 = (beta^(t-1))*( pxest(2,m)^(t-1));
temp3 = (beta^(t-1))*( (pxest(2,m)^(t-2))*(t-1));
B = B+temp;
v_crime_bar(m) = v_crime_bar(m) + temp*extp(m);
gamma(m) = gamma(m) + temp2;
end
end
v_crime_bar = v_crime_bar/B;
v_crime_bar_big(:,val) = v_crime_bar;
end

gammaar2 = v_crime_bar_big(:,2) - v_crime_bar_big(:,1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load est_res_s2 rentd  M con_m
load datafromstata gamma 

statmwtp = rentd;

picmwtpar2 = 0*rentd;
for m = 1:M
    index = find(con_m==m);
    picmwtpar2(index) = statmwtp(index)/gammaar2(m);
end

for m = 1:M
    index = find(con_m==m);
    res2(m,:) = [ mean(picmwtpar2(index))];
end

rob_AR2 = res2;
save rob_AR2 rob_AR2

